clc
clear
% close all
dbstop if error
% 设定控制条件
residual_newton_max = 0.0001;
residual_newton = 1E5;
newton_count_max = 20;
% 设定初始条件
x = -1:0.001:1;  % 区域,间断在0
t = 0:0.01:0.14;
z = 0:0.001:0.4;
global gamma p1 rho1 c1 p2 rho2 c2 u1 u2;

gamma = 1.2;
p1 = 1;
rho1 = 1;
c1 = sqrt(gamma*p1/rho1);
p2 = 0.1;
rho2 = 0.125;
c2 = sqrt(gamma*p2/rho2);
u1 = 0;
u2 = 0;
R = 8.314;

p_star = 0:0.001:2;
y = zeros(1,length(p_star));
for i = 1:length(p_star)

    y(i) = F(p_star(i)) + u1 - u2;
end
figure(1);
plot(p_star,y,'LineWidth',2);
legend("F")
y = zeros(length(x),length(t),4);
% newton迭代求解p_star;
% 首先给定初值
p_star0 = (p1+p2)/2;
p_star = p_star0;
delta_p_star = 0;
newton_count = 0;
while(residual_newton > residual_newton_max && newton_count < newton_count_max)
    delta_p_star = -F(p_star)/F_diff(p_star);
    p_star = p_star + delta_p_star;
    if (newton_count == 0)
        p_star0 = p_star;
    end
    newton_count=newton_count+1 ;
    residual_newton = abs(delta_p_star)/p_star0;
end
% 得到p_star,求出star区域的其他物理量
if(p_star > p1)
    % 两侧均为激波
end
if(p_star>p2 && p_star < p1)
    %左膨胀波，右侧激波
    u_star = u1 - f1(p_star);
    rho_star_L = (p_star*rho1^gamma/p1)^(1/gamma);
    cL = sqrt(gamma*p_star/rho_star_L);
    Z2 = u2 -(p_star -p2)/(rho2*u2 - rho2*u_star);
    rho_star_R = rho2*(u2 -Z2)/(u_star - Z2);
    cR = sqrt(gamma*p_star/rho_star_R);
    % 还有左侧稀疏波中的物理量
    % 确定各区的范围
    for j =1:length(t)
        for i =1:length(x)
            % 在左侧未扰动区
            if(x(i) < (u1-c1)*t(j))
                y(i,j,1) = u1;
                y(i,j,2) = p1;
                y(i,j,3) = rho1;
                y(i,j,4) = c1^2/gamma/R;
            end
            % 在稀疏波内
            if(x(i) < (u_star -cL)*t(j) && x(i) >= (u1-c1)*t(j))
                %不同位置波速是不一样的，稀疏波有厚度
                c_expansion = (gamma-1)/(gamma+1)*(u1-x(i)/t(j))+2/(gamma+1)*c1;
                u_expansion = c_expansion + x(i)/t(j);
                p_expansion = p1*(c_expansion/c1)^(2*gamma/(gamma-1));
                rho_expansion = gamma*p_expansion/c_expansion^2;
                y(i,j,1) = u_expansion;
                y(i,j,2) = p_expansion;
                y(i,j,3) = rho_expansion;
                y(i,j,3) = c_expansion^2/gamma/R;
            end
            % 在中间 接触间断左侧
            if(x(i) <u_star*t(j) && x(i) >= (u_star-cL)*t(j))
                y(i,j,1) = u_star;
                y(i,j,2) = p_star;
                y(i,j,3) = rho_star_L;
                y(i,j,4) = cL^2/gamma/R;
            end


            % 在中间 接触间断右侧
            if(x(i) >=u_star*t(j) && x(i) < Z2*t(j))
                y(i,j,1) = u_star;
                y(i,j,2) = p_star;
                y(i,j,3) = rho_star_R;
                y(i,j,4) = cR^2/gamma/R;
            end

            % 在右侧未扰动区 激波波前
            if(x(i) >= Z2*t(j))
                y(i,j,1) = u2;
                y(i,j,2) = p2;
                y(i,j,3) = rho2;
                y(i,j,4) = c2^2/gamma/R;
            end
        end
    end
end
if(p2>p_star)
    %两侧膨胀波
end
if(p_star<=0)
    % 中间真空
end

% 展示结果

figure(2)
test_data = load("lixinliang.txt");
plot(x, y(:,end,2), 'LineWidth',1);
hold on
plot(test_data(:,1) , test_data(:,2), 'o', 'LineWidth',1)
legend("Present", "Li");
xlabel("x");
ylabel("p");
axis([-1 1 0 1.1])



figure(3);
% set(gcf,'outerposition',get(0,'screensize'));
[X,Z] = meshgrid(x,z);

[row, col] = size(X);
result2D = zeros(row, col,length(t), 4);
for j =1:length(t)
    for i =1:row
        for k = 1:col
            % 在左侧未扰动区
            if(X(i,k) < (u1-c1)*t(j))
                result2D(i,k,j,1) = u1;
                result2D(i,k,j,2) = p1;
                result2D(i,k,j,3) = rho1;
                result2D(i,k,j,4) = c1^2/gamma/R;
            end
            % 在稀疏波内
            if(X(i,k) < (u_star -cL)*t(j) && X(i,k) >= (u1-c1)*t(j))
                %不同位置波速是不一样的，稀疏波有厚度
                c_expansion = (gamma-1)/(gamma+1)*(u1-X(i,k)/t(j))+2/(gamma+1)*c1;
                u_expansion = c_expansion + X(i,k)/t(j);
                p_expansion = p1*(c_expansion/c1)^(2*gamma/(gamma-1));
                rho_expansion = gamma*p_expansion/c_expansion^2;
                result2D(i,k,j,1) = u_expansion;
                result2D(i,k,j,2) = p_expansion;
                result2D(i,k,j,3) = rho_expansion;
                result2D(i,k,j,4) = c_expansion^2/gamma/R;
            end
            % 在中间 接触间断左侧
            if(X(i,k) <u_star*t(j) && X(i,k) >= (u_star-cL)*t(j))
                result2D(i,k,j,1) = u_star;
                result2D(i,k,j,2) = p_star;
                result2D(i,k,j,3) = rho_star_L;
                result2D(i,k,j,4) = cL^2/gamma/R;
            end


            % 在中间 接触间断右侧
            if(X(i,k) >=u_star*t(j) && X(i,k) < Z2*t(j))
                result2D(i,k,j,1) = u_star;
                result2D(i,k,j,2) = p_star;
                result2D(i,k,j,3) = rho_star_R;
                result2D(i,k,j,4) = cR^2/gamma/R;
            end

            % 在右侧未扰动区 激波波前
            if(X(i,k) >= Z2*t(j))
                result2D(i,k,j,1) = u2;
                result2D(i,k,j,2) = p2;
                result2D(i,k,j,3) = rho2;
                result2D(i,k,j,4) = c2^2/gamma/R;
            end
        end
    end
end

makeVideo = VideoWriter("ShockTube","MPEG-4");
makeVideo.FrameRate = 5;
open(makeVideo);

for i=1:length(t)

    [p,s] = contourf(X,Z,result2D(:,:,i,3),100);
    xlabel("x");
    ylabel("r");

    shading interp;
    s.LineStyle = 'none';
    axis equal;
    colormap  jet;
    colorbar;
    drawnow;
    frame = getframe(gcf);
    writeVideo(makeVideo,frame);
    ylabel(colorbar,'\rho','FontName','Times New Roman','FontWeight','bold')
    %     pause(0.2)
end
close(makeVideo);

